module compute
    use basicmodule
	use setthmod
	use analyticalIS
	use credmod
    implicit none
	
	integer				:: nsimeffmax
	
	contains
	

	subroutine prep
		real, allocatable	:: eps(:,:,:)
		real	:: Y(0:k),h,th(nth_one)
		integer	:: ims(npb,nblocks)
		integer	:: l,lb,i,ic

		allocate(eps(0:k,npb,nblocks))
		
		do lb=1,nblocks
			call rnund(nimp,ims(:,lb))
		enddo
		do lb=1,nblocks
			call rnnoa(eps(0,:,lb))
			do l=1,npb
				call rnun(eps(1:k,l,lb))
				eps(1:k,l,lb)=-log(eps(1:k,l,lb))
			enddo
		enddo
!$omp parallel do private(l,th,Y,h,i)		
		do lb=1,nblocks
			do l=1,npb
				th=thlistimp(:,ims(l,lb))
				call setYfromeps(th,eps(:,l,lb),Y)
				Y0s(:,l,lb)=Y
				h=0
				do i=1,nimp
					h=h+getdensimp(thlistimp(:,i),Y)
				enddo
				densimp(l,lb)=h/nimp
				falt(l,lb)=getf1(Y(1:))
			enddo
		enddo
		print *,"done with prep"
		call printtime
	end subroutine
	
	subroutine checkIS
		real	:: LR(4,nimp,2),th(3),y(0:k),c(4)
		integer	:: lb,l,i
		real	:: x(3),thlist(3,nimp)
		
		do i=1,nimp
			call rnun(x)
			thlist(:,i)=get_msx0(x)
		enddo
		thlist=thlistimp
		LR=0
!$omp parallel do private(l,i,y,c) reduction(+:LR)
		do lb=1,nblocks
			do l=1,npb
				y=Y0s(:,l,lb)
				do i=1,nimp
					c=getdens0(thlist(:,i),y)/densimp(l,lb)*[1,boole(y(1)<0),boole(y(1)/y(k)<swflat)*boole(y(k)>0),boole(y(1)/y(k)<swflat)*boole(y(k)>0)*boole(y(1)>swsmall)]
					if(.not. all(isfinite(c))) then
						c=getdens0(thlist(:,i),y)
					endif
					LR(:,i,1)=LR(:,i,1)+c
					LR(:,i,2)=LR(:,i,2)+c**2
					
				enddo
			enddo
		enddo
		LR=LR/nsim
		call mdisp(LR(:,:,1))
		LR(:,:,2)=sqrt((LR(:,:,2)-LR(:,:,1)**2)/nsim)
		call mdisp(LR(:,:,2))
		print *,"IS=1 t-stats"
		call mdisp((LR(1,:,1)-1)/LR(1,:,2))
			
	end subroutine
	
	subroutine preppairs
		integer	:: lb1,l1,l2,lb2,i2,ii
		integer	:: ix
		integer(8)	:: code
		
		nsimeffmax=0
		if(allocated(pairsdat)) deallocate(pairsdat,pairsind)
		allocate(pairsdat(nsimeffmax),pairsind(4,nsimeffmax))
		nsimeff=0
		do lb1=1,nblocks
			lb2=lb1+1
			if(lb2>nblocks) lb2=1
			pairbind(lb1)=nsimeff
			do l1=1,npb
				l2=0
				do i2=1,nbi
					code=credflags(i2,l1,lb1)
					do ii=1,nfpc
						l2=l2+1
						if(code<0) then
							call addx(l1,lb1,l2,lb2)
						endif
						code=2*code
					enddo
				enddo
			enddo
		enddo
		print *,"done with preppairs"
		print *,"nsimeff is",nsimeff
		call printtime
		call modifypairs
		call printtime
	end subroutine
	
	subroutine modifypairs
		integer	:: lb,l1,i2,ix,l,j
		integer(8)	:: ccode,tcode
		integer	:: rejind
		logical	:: rej
		
		if(.not. tstconst) then
			ctest=credflags
			return
		endif

!$omp parallel do private(l,l1,i2,ix,ccode,tcode,rejind,rej)
		do lb=1,nblocks
			l=pairbind(lb)
			do l1=1,npb
				do i2=1,nbi
					ccode=credflags(i2,l1,lb)
					tcode=0
					do ix=1,nfpc
						tcode=tcode*2
						if(ccode<0) then
							l=l+1
							tcode=tcode+1
							if(level/=0.05) then
								rejind=0
								rej=rejectfrompair(pairsdat(l),tst5)
								if(level>0.05 .and. rej) then
									rejind=1
								elseif(level<0.05 .and. (.not. rej)) then
									rejind=-1
								endif
								if(rejind==0) then
									if(level>0.1) then
										if(rejectfrompair(pairsdat(l),tst10)) rejind=1
									elseif(level<0.05 .and. level>0.01) then
										if(rejectfrompair(pairsdat(l),tst1)) rejind=1
									endif
								endif
								if(rejind>0) then
									pairsdat(l)%f1s=1E200; pairsdat(l)%f1b=1E200
								elseif(rejind<0) then
									pairsdat(l)%f1s=1E-200; pairsdat(l)%f1b=1E-200
									tcode=tcode-1
								endif
							endif
						endif
						ccode=2*ccode
					enddo
					ctest(i2,l1,lb)=tcode*2**(64-nfpc)
				enddo
			enddo
		enddo
	end subroutine
	
	subroutine addx(l1,lb1,l2,lb2)
		integer	:: l1,lb1,l2,lb2
		type(ptype), allocatable	:: pairsdatnew(:)
		type(ptype)			:: p
		integer,allocatable		:: pairsindnew(:,:)
		logical	:: rej
		integer			:: rejind
			
		if(nsimeff==nsimeffmax) then
			nsimeffmax=nsimeff+5000000
			allocate(pairsdatnew(nsimeffmax),pairsindnew(4,nsimeffmax))
			pairsdatnew(1:nsimeff)=pairsdat
			pairsindnew(:,1:nsimeff)=pairsind
			call move_alloc(pairsdatnew,pairsdat)
			call move_alloc(pairsindnew,pairsind)
		endif
		call setpair(Y0s(:,l1,lb1),Y0s(:,l2,lb2),p)
		p%f1s=falt(l1,lb1)
		p%f1b=falt(l1,lb1)*falt(l2,lb2)

		nsimeff=nsimeff+1
		pairsdat(nsimeff)=p
		pairsind(:,nsimeff)=[l1,lb1,l2,lb2]
	end subroutine
	
	subroutine setpair(y1,y2,p)
		real			:: y1(0:k),y2(0:k)
		type(ptype)		:: p
		p%Z(:,1)=y1(1:k)
		p%Z(:,2)=y2(1:k)
		p%Zsum=y1(0)-y2(0)
	end subroutine
	
	function filename() result(val)
		character(len=50)	:: val
		val=convtos(k)//"_"//convtos(n0)//"_"//convtos(nint(1000*level))//".txt"
	end function
	
	subroutine savelam
		if(stage==1) then
			lam_stage1=lam
			thlist0_stage1=thlist0(1:3,1:n0active)
			n0active_stage1=n0active
			call savemat(trim(folder)//"lamth_stage1_"//trim(filename()),lam_stage1.cud.thlist0_stage1)
		else			
			call savemat(trim(folder)//"lamth_"//trim(filename()),lam(1:n0active).cud.thlist0(:,1:n0active))
		endif
	end subroutine
	
	subroutine loadlam
		real, allocatable	:: dta(:,:)
		dta=loadmat(trim(folder)//"lamth_stage1_"//trim(filename()))
		lam_stage1=dta(1,:)
		thlist0_stage1=dta(2:4,:)
		n0active_stage1=size(lam_stage1)

		dta=loadmat(trim(folder)//"lamth_"//trim(filename()))
		n0active=size(dta,2)
		lam=dta(1,:)
		thlist0(:,1:n0active)=dta(2:7,:)
		call mdisp(lam.cud.thlist0(:,1:n0active))
		call setcvs
	end subroutine
	
	subroutine setfintest
		integer	:: lb,l1,i2,ix,l,j
		integer(8)	:: ccode,tcode
		real		:: h
!$omp parallel do private(l,l1,i2,ix,tcode,ccode,j,h)
		do lb=1,nblocks
			l=pairbind(lb)
			do l1=1,npb
				do i2=1,nbi
					tcode=0
					ccode=credflags(i2,l1,lb)
					do ix=1,nfpc
						tcode=2*tcode
						if(ccode<0) then
							l=l+1
							h=0
							do j=1,n0active
								h=h+lam(j)*getdens0p(thlist0(:,j),pairsdat(l))
							enddo
							if(h<cvglobal) tcode=tcode+1
						endif
						ccode=2*ccode
					enddo
					ctest(i2,l1,lb)=tcode*2**(64-nfpc)
				enddo
			enddo
		enddo
	end subroutine
	
	function rejectfrompair(p,tst) result(val)
		type(ptype)	:: p
		type(tsttype)	:: tst
		logical	:: val
		real	:: y1(0:k),y2(0:k),fa1,fa2
		y1(0)=p%Zsum
		y1(1:k)=p%Z(:,1)
		y2(0)=0
		y2(1:k)=p%Z(:,2)
		fa1=p%f1s
		fa2=p%f1b/fa1
		val=rejectfast(y1,y2,fa1,fa2,tst)
	end function
		
	function rejectfast(y1,y2,fa1,fa2,tst) result(val)
		type(tsttype)	:: tst
		real	:: y1(0:k),y2(0:k),fa1,fa2
		logical	:: val
		type(ptype)	:: p
		real	:: h
		integer	:: i
		
		if(.not. getapriorirej(y1,y2,fa1,fa2,tst)) then
			val=.false.
			return
		endif

		if(tst%stage==1) then
			val=rej_stage1(y1,y2,fa1,tst)
			return
		endif

		call setpair(Y1,Y2,p)
		p%f1s=fa1
		p%f1b=fa1*fa2 

		h=0
		do i=1,tst%n0active
			h=h+tst%lam(i)*getdens0p_stage2(tst%thlist0(:,i),p)
		enddo
		h=h/p%f1b
		val=h<cvglobal
	end function

	function reject(Yo,tst) result(val)
		type(tsttype)	:: tst
		real	:: Yo(0:2*k)
		logical	:: val
		real	:: Ys(0:k,2),y1b(2),y(0:k),fas(2)
		logical	:: rej
		
		!if(Yo(1)<0 .or. Yo(k+1)<0) then
		!	val=.true.
		!	return
		!endif
		Ys(:,1)=Yo(0:k)
		Ys(:,2)=[0.0,Yo(k+1:)]
		fas=[getf1(Yo(1:k)),getf1(Yo(k+1:))]
		if(tstconst) then
			if(level/=0.05) then
				rej=rejectfast(Ys(:,1),Ys(:,2),fas(1),fas(2),tst5)
				if(level>0.05 .and. rej) then
					val=.true.
					return
				elseif(level<0.05 .and. (.not. rej)) then
					val=.false.
					return
				endif
				if(level>0.1) then
					if(rejectfast(Ys(:,1),Ys(:,2),fas(1),fas(2),tst10)) then
						val=.true.; return
					endif
				elseif(level<0.05 .and. level>0.01) then
					if(rejectfast(Ys(:,1),Ys(:,2),fas(1),fas(2),tst1)) then
						val=.true.; return
					endif
				endif
			endif
		endif
		val=rejectfast(Ys(:,1),Ys(:,2),fas(1),fas(2),tst)
	end function		
	
end module
